---
title: "Analysis of the Skill-Based Test"
output:
  pdf_document: default
  html_document: default
date: "2025-10-06"
---

## Packages used:

```{r, echo = T, results = 'hide', error=FALSE, warning=FALSE, message=FALSE}
rm(list = ls())
library(tidyverse)
library(ggplot2)
library(ggrepel)
library(ggpubr)
library(estimatr)
library(dotwhisker)
library(broom)
library(dplyr)
library(sjPlot)
library(hrbrthemes)
library(ivreg)
library(readxl)
library(gtools)
library(stargazer)
library(cobalt)
library(summarytools)
library(pwrss)
library(pwr)
```

## Importing, Cleaning, & Stacking Data. 

```{r,echo = T, results = 'hide',  error=FALSE, warning=FALSE, message=FALSE}
milexp1 <- read_xlsx("milexp1.xlsx")
milexp2 <- read_xlsx("milexp2.xlsx")
```

```{r,echo = T, results = 'hide',  error=FALSE, warning=FALSE, message=FALSE}
milexp1$T1 = rep(NA, nrow(milexp1))
milexp1$T1[milexp1$group == 1] <- 1 # partisan agreement
milexp1$T1[milexp1$group == 2] <- 1 # partisan agreement
milexp1$T1[milexp1$group == 3] <- 0 # partisan disagreement
milexp1$T1[milexp1$group == 4] <- 0 # partisan disagreement

milexp1$T2 = rep(NA, nrow(milexp1))
milexp1$T2[milexp1$group == 1] <- 1 # incumbent party identifier
milexp1$T2[milexp1$group == 2] <- 0 # opposition party identifier
milexp1$T2[milexp1$group == 3] <- 1 # incumbent party identifier
milexp1$T2[milexp1$group == 4] <- 0 # opposition party identifier

milexp1$female[milexp1$female == 1] <- 2
milexp1$female[milexp1$female == 0] <- 1

milexp1$priorexp[milexp1$priorexp == 0] <- 2

milexp1$study = rep(NA, nrow(milexp1))
milexp1$study <- 1
milexp1$ID <- (1:nrow(milexp1))

milexp1 <- rename(milexp1, "risklife" = "SLE1-risklife")
milexp1 <- rename(milexp1, "risksubordin" = "SLE2-risksubordin")
milexp1 <- rename(milexp1, "discretion" = "SLE3-discretion")
milexp1 <- rename(milexp1, "interunit" = "SLE4-interunit")
milexp1 <- rename(milexp1, "unitcohesion" = "SLE5-unitcohesion")
milexp1 <- rename(milexp1, "obey" = "SLE6-obey")

milexp1 <- rename(milexp1, "AM_1" = "AM1-threat")
milexp1 <- rename(milexp1, "AM_2" = "AM2-partytype")
milexp1 <- rename(milexp1, "AM_3" = "AM3-artificial")


milexp2$T1 = rep(NA, nrow(milexp2))
milexp2$T1[milexp2$group == 1] <- 1 # partisan agreement
milexp2$T1[milexp2$group == 2] <- 1 # partisan agreement
milexp2$T1[milexp2$group == 3] <- 0 # partisan disagreement
milexp2$T1[milexp2$group == 4] <- 0 # partisan disagreement

milexp2$T2 = rep(NA, nrow(milexp2))
milexp2$T2[milexp2$group == 1] <- 1 # incumbent party identifier
milexp2$T2[milexp2$group == 2] <- 0 # opposition party identifier
milexp2$T2[milexp2$group == 3] <- 1 # incumbent party identifier
milexp2$T2[milexp2$group == 4] <- 0 # opposition party identifier

milexp2$female[milexp2$female == 1] <- 2
milexp2$female[milexp2$female == 0] <- 1

milexp2$study = rep(NA, nrow(milexp2))
milexp2$study <- 2
milexp2$ID <- ((nrow(milexp1)+1):(nrow(milexp2)+nrow(milexp1)))

milexp2 <- rename(milexp2, "risklife" = "SLE1-risklife")
milexp2 <- rename(milexp2, "risksubordin" = "SLE2-risksubordin")
milexp2 <- rename(milexp2, "discretion" = "SLE3-discretion")
milexp2 <- rename(milexp2, "interunit" = "SLE4-interunit")
milexp2 <- rename(milexp2, "unitcohesion" = "SLE5-unitcohesion")
milexp2 <- rename(milexp2, "obey" = "SLE6-obey")

milexp2 <- rename(milexp2, "AM_1" = "AM1-threat")
milexp2 <- rename(milexp2, "AM_2" = "AM2-partytype")
milexp2 <- rename(milexp2, "AM_3" = "AM3-artificial")
```

```{r, error=FALSE, warning=FALSE, message=FALSE}
milexp1.1 <- milexp1 %>% 
  select(risklife:obey, NIM:partyIDstrength, anger:ID, group) %>%
  group_by(ID)
milexp2.1 <- milexp2 %>% 
  select(NIM:AM_3, age:partyIDstrength, anger:group, T1:ID) %>%
  group_by(ID)

milexp <- rbind(milexp1.1, milexp2.1)

datstack <- milexp %>%
  pivot_longer(cols = "risklife":"obey", 
               names_to = "Variable", 
               values_to = "Value",
               values_drop_na = F)
```

Testing H1 with behavioral measures (shooting scores): We report the effect of partisan disagreement on shooting scores. Although the sample size is too small, the analysis is underpowered, and (maybe due to that) the correlation is not statistically significant, partisan agreement is positively correlated with shooting score. We provide results from simple OLS regression (with plots included), permutation test, and ex post power analysis.

```{r, error=FALSE, warning=FALSE, message=FALSE}
set.seed(123)

## Difference in Means

milexp1.2 <- milexp1[complete.cases(milexp1$shooting), ]


mean(milexp1.2$shooting)
sd(milexp1.2$shooting)

mean(milexp1.2$shooting[milexp1.2$T1 == "1"])
sd(milexp1.2$shooting[milexp1.2$T1 == "1"])

mean(milexp1.2$shooting[milexp1.2$T1 == "0"])
sd(milexp1.2$shooting[milexp1.2$T1 == "0"])


test.stat1 <- abs(mean(milexp1.2$shooting[milexp1.2$T1 == "1"]) - 
                  mean(milexp1.2$shooting[milexp1.2$T1 == "0"])) 

test.stat1


## OLS (Shooting Score as DV)

H13 <- lm(shooting ~ T1, data = milexp1)

summary(H13)

H14 <- lm(shooting ~ T1 + age + female + hometown + anger + risk + moral + dove,
          data = milexp1)

summary(H14)


## The effect of partisan disagreement

milexp1$Elite[milexp1$T1 == 0] <- "Disagree"
milexp1$Elite[milexp1$T1 == 1] <- "Agree"
milexp1$Partisanship[milexp1$T2 == 0] <- "Opposition"
milexp1$Partisanship[milexp1$T2 == 1] <- "Incumbent"

milexp1 %>%
  mutate(name = fct_relevel(Elite, 
            "Agree", "Disagree")) %>%
  ggplot(aes(x = shooting, fill = name)) +
  geom_density(aes(fill = name), alpha = .3) + theme_bw() + xlab("Shooting Score") +
  ylab("Density") + geom_vline(xintercept = mean(milexp1$shooting[milexp1$T1 == 0], na.rm = TRUE),
                               colour = "pink3", linetype = 2) +
  geom_vline(xintercept = mean(milexp1$shooting[milexp1$T1 == 1], na.rm = TRUE),
                               colour = "skyblue3", linetype = 2) +
  scale_fill_manual(name = "Elite", values=c("skyblue3", "pink3"))


## The effect of partisan disagreement among Opposition party members
milexp1 %>%
  filter(T2 == 0) %>%
  mutate(name = fct_relevel(Elite, 
            "Agree", "Disagree")) %>%
  ggplot(aes(x = shooting, fill = name)) +
  geom_density(aes(fill = name), alpha = .3) + theme_bw() +
  xlab("Shooting Score (among the opposition party identifiers)") +
  ylab("Density") + geom_vline(xintercept = mean(milexp1$shooting[milexp1$T1 == 1 & milexp1$T2 == 0],
                                                 na.rm = TRUE), colour = "skyblue3", linetype = 2) +
  geom_vline(xintercept = mean(milexp1$shooting[milexp1$T1 == 0 & milexp$T2 == 0],
                               na.rm = TRUE), colour = "pink3", linetype = 2) +
  scale_fill_manual(name = "Elite", values=c("skyblue3", "pink3"))


## Permutation test

n <- length(milexp1.2$shooting)
P <- 100000
PermSamples <- matrix(0, nrow = n, ncol = P)
shooting <- milexp1.2$shooting

for(i in 1:P)
  {
    PermSamples[, i] <- sample(shooting, 
                               size = n, 
                               replace = FALSE)
  }

Perm.test.stat1 <- rep(0, P)

for (i in 1:P)
  {
    Perm.test.stat1[i] <- abs(mean(PermSamples[milexp1.2$T1 == "1",i]) - 
                              mean(PermSamples[milexp1.2$T1 == "0",i]))
  }

mean(round(Perm.test.stat1, 2))

# Permutation p-value

mean(Perm.test.stat1 >= test.stat1)



## Ex Post Power Analysis
# Power of H13 analysis: About .24
power.t.test(ncp = 1.24, df = 277, alpha = .05, alternative = "two.sided",
             plot = T)

# n > 1310 necessary to achieve the statistical power of .8
# for detecting the effect at the .95 significance level,
# according to the current finding
pwr.f2.test(u = 1, sig.level = .05, power = .8, f2 = .006)
# f2 based on the R-squared value from H13
``` 